function triangular_systems()
% TRIANGULAR_SYSTEMS 三角系统求解演示
% 
% 演示上三角和下三角系统的前向和后向替换算法

fprintf('=== 三角系统求解 ===\n\n');

%% 1. 下三角系统 - 前向替换
fprintf('1. 下三角系统求解 (前向替换)\n');
fprintf('求解 L*x = b，其中L是下三角矩阵\n\n');

% 创建下三角矩阵
L = [2, 0, 0; 
     1, 3, 0; 
     4, 2, 5];
b = [4; 7; 20];

fprintf('下三角矩阵 L =\n');
disp(L);
fprintf('右端向量 b = [%g; %g; %g]\n', b);

% 前向替换算法
x_forward = forward_substitution(L, b);
fprintf('前向替换解: x = [%.6f; %.6f; %.6f]\n', x_forward);

% 验证解
residual = L * x_forward - b;
fprintf('残差 ||L*x - b|| = %e\n', norm(residual));

% 与MATLAB内置求解器比较
x_matlab = L \ b;
fprintf('MATLAB解: x = [%.6f; %.6f; %.6f]\n', x_matlab);
fprintf('差异: ||x_forward - x_matlab|| = %e\n', norm(x_forward - x_matlab));

%% 2. 上三角系统 - 后向替换
fprintf('\n2. 上三角系统求解 (后向替换)\n');
fprintf('求解 U*x = b，其中U是上三角矩阵\n\n');

% 创建上三角矩阵
U = [3, 2, 1; 
     0, 4, 3; 
     0, 0, 2];
b2 = [14; 20; 6];

fprintf('上三角矩阵 U =\n');
disp(U);
fprintf('右端向量 b = [%g; %g; %g]\n', b2);

% 后向替换算法
x_backward = backward_substitution(U, b2);
fprintf('后向替换解: x = [%.6f; %.6f; %.6f]\n', x_backward);

% 验证解
residual2 = U * x_backward - b2;
fprintf('残差 ||U*x - b|| = %e\n', norm(residual2));

% 与MATLAB内置求解器比较
x_matlab2 = U \ b2;
fprintf('MATLAB解: x = [%.6f; %.6f; %.6f]\n', x_matlab2);
fprintf('差异: ||x_backward - x_matlab|| = %e\n', norm(x_backward - x_matlab2));

%% 3. 算法复杂度分析
fprintf('\n3. 算法复杂度分析\n');
fprintf('三角系统求解的时间复杂度为 O(n^2)\n');

% 测试不同规模矩阵的求解时间
sizes = [100, 200, 500, 1000];
times_forward = zeros(size(sizes));
times_backward = zeros(size(sizes));

fprintf('\n矩阵规模    前向替换时间    后向替换时间\n');
fprintf('--------    ------------    ------------\n');

for i = 1:length(sizes)
    n = sizes(i);
    
    % 生成随机三角矩阵
    L_test = tril(randn(n, n)) + eye(n);  % 确保对角元非零
    U_test = triu(randn(n, n)) + eye(n);
    b_test = randn(n, 1);
    
    % 测试前向替换
    tic;
    x_test1 = forward_substitution(L_test, b_test);
    times_forward(i) = toc;
    
    % 测试后向替换
    tic;
    x_test2 = backward_substitution(U_test, b_test);
    times_backward(i) = toc;
    
    fprintf('%6d      %8.4f ms      %8.4f ms\n', n, times_forward(i)*1000, times_backward(i)*1000);
end

%% 4. 病态三角系统
fprintf('\n4. 病态三角系统\n');
fprintf('当三角矩阵的对角元素很小时，系统变得病态\n');

% 创建病态下三角矩阵
L_ill = [1e-10, 0, 0; 
         1, 1, 0; 
         1, 1, 1];
b_ill = [1e-10; 2; 3];

fprintf('病态下三角矩阵 L =\n');
disp(L_ill);
fprintf('条件数: %.2e\n', cond(L_ill));

x_ill = forward_substitution(L_ill, b_ill);
fprintf('解: x = [%.6e; %.6f; %.6f]\n', x_ill);

% 添加小扰动
delta_b = [1e-12; 0; 0];
b_perturbed = b_ill + delta_b;
x_perturbed = forward_substitution(L_ill, b_perturbed);

rel_error_b = norm(delta_b) / norm(b_ill);
rel_error_x = norm(x_perturbed - x_ill) / norm(x_ill);

fprintf('右端扰动的相对误差: %e\n', rel_error_b);
fprintf('解的相对误差: %e\n', rel_error_x);
fprintf('误差放大因子: %.2e\n', rel_error_x / rel_error_b);

end

function x = forward_substitution(L, b)
% FORWARD_SUBSTITUTION 前向替换算法
% 求解下三角系统 L*x = b
%
% 输入:
%   L - n×n 下三角矩阵
%   b - n×1 右端向量
% 输出:
%   x - n×1 解向量

n = size(L, 1);
x = zeros(n, 1);

for i = 1:n
    sum_term = 0;
    for j = 1:i-1
        sum_term = sum_term + L(i, j) * x(j);
    end
    x(i) = (b(i) - sum_term) / L(i, i);
end

end

function x = backward_substitution(U, b)
% BACKWARD_SUBSTITUTION 后向替换算法
% 求解上三角系统 U*x = b
%
% 输入:
%   U - n×n 上三角矩阵
%   b - n×1 右端向量
% 输出:
%   x - n×1 解向量

n = size(U, 1);
x = zeros(n, 1);

for i = n:-1:1
    sum_term = 0;
    for j = i+1:n
        sum_term = sum_term + U(i, j) * x(j);
    end
    x(i) = (b(i) - sum_term) / U(i, i);
end

end